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ABSTRACT 

We explore some properties of twin kilohertz quasiperiodic oscillations (QPOs) in a simple toy-model consisting of two oscillation 
modes coupled by a general nonlinear force. We examine resonant effects by slowly varying the values of the tunable, and nearly 
commensurable, eigenfrequencies. The behavior of the actual oscillation frequencies and amplitudes during a slow transition through 
the 3:2 resonance is examined in detail and it is shown that both are significantly affected by the nonlinearities in the governing 
equations. In particular, the amplitudes of oscillations reflect a resonant exchange of energy between the modes, as a result the 
initially weaker mode may become dominant after the transition. We note that a qualitatively similar behavior has been recently 
reported in several neutron star sources by Torok (2008), who found that the difference of amplitudes in neutron star twin peak QPOs 
changes sign as the observed frequency ratio of the QPOs passes through the value 3:2. 
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> 1. Introduction 

, There is little doubt that kHz QPOs observed in low mass X- 
C^) • ray binaries originate in the inner region of accretion disks of 
CO ' the neutron star or black hole sources. The very fact that these 
' \ quasi-periodic modulations are observed in the X-ray flux, i.e., 
f~>) ■ in the most luminous part of the source spectrum — and at the 
' same time at high photon energies — implies an origin in the in- 
f^i , ner part of the accretion flow, where most of the gravitational po- 
I* ■ tential energy of the accreting matter is released. Time-of-flight 
. £^ | arguments also imply that a kHz modulation of unbeamed radi- 
, ation must originate in a region not larger than ~ 100 km, with 
5_i ' an even sharper bound when the high degree of coherence of the 
observed modulation is taken into account. 

In the present paper we focus on kHz QPOs that show u p in 
pairs, the so called twin peak QPOs (see Ivan der Klis 2000, for 
a review), and we assume that they correspond to oscillations of 
an accretion-disk structure. Reproducible features of these high- 
frequency QPOs imply that they reflect some fundamental prop- 
erties of the systems, such as their mass and spin. 

Many mo dels have been proposed, some involvin g or- 
bital motions (IKluzniak et al.lll990t Istella & Vietril 1 1 9991)7 oth- 
ers based on oscillations of the accretion disk (I Wagoner! 



Petri! 12005b iPetril 120061) and of tori dRezzolla et al.l l2003blat 



19991 Wag oner et al. 1]200U !Katoll200T 



Bursa et al.ll2004b iMontero et alj|2004 



2004al l2005afbl 12008 



(V 



2007b iBlaes et al J 120061 



2007b ISchnittman & Rezzollall2006b IZanotti et al.ll2005b faora 



2008). Most of these models involve two or more oscillatory mo- 



tions, which may interact in a non-linear manner. 



The appearance of rational ratios (mainly ~3:2) be- 
tween the centroid of twin peak QPOs points at a non- 
linear coupling between diffe rent modes of oscillations 
(lAbramowicz & Kluzniakl 1200 lb . The imp ortance of a res- 
onant p rocess was also pointed out b y iPsaltis & Norman! 
(12000b: IKluzniak & Abramowiczl (|2001|) : iTitarchukl d2002l) : 
IKluzniak et al.l (l2004f) and in v estigated i n mor e detail, e.g., 

iHorak et al] 



in Abramowicz et al. 



(l2004l) : lLee et al.l (12004): Hor ak & Karasl (|200< 



2003bl): iRebuscol 
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One example of an investigation of resonances between spe- 
cific modes in an accretion disk can be found in the work o f lKatol 
(l2005bl) . who proposed a model in which the two kHz QPOs are 
due to a nonlinear resonant interaction between the oscillations 
of a thin accretion disk and an oscillatory perturbation induced 
by the non-axisymmetric rotating neutron star. Both kHz oscil- 
lations are localized close to a resonant radius, where the inter- 
action with the perturbation is most prominent. The spatial prox- 
imity of the two oscillations may facilitate an internal resonan ce 
between them. As pointed out recently by iMao et ail d2008l) . a 
resonance may occur between two modes with non-overlapping 
wavefunctions, when it is mediated by a traveling wave. In other 
words, the two oscillators in resonance need not be located in the 
same annulus of the accretion disk. More recently iKatol (120081) 
relaxed the local approximation and analyzed the nonlinear reso- 
nant excitation of inertial-acoustic disk oscillations in deformed 
disks (a warped disk or a one-armed pattern symmetric with re- 
spect to the equator). 

Other QPO models involv e vortices in the d i sk or mag- 
netic interactions. For example iTagger & Varniera (|2006j) sug- 
gest the formation of t wo or three Rossby-w ave vortices in mi- 
croquasar disks, while iLi & Naravanl d2004l) suggested a model 
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in which the two QPOs are identified with two Rayleigh-Taylor 
or Kelvin-Helmholtz unstable modes developed at the inter- 
face between the disk and the magnetosphere. Similarly, in the 
Alfven wave oscillations model dZhangl2004HZhang et al.l2007l) 
the two QPOs correspond to the Keplerian orbital frequency and 
the frequency of the Alfven waves excited at the same radius. 
It is less clear whether the resonant interaction described in this 
paper would apply to those models as well. 

In this paper we consider some general properties of non- 
linear resonance that may be applicable to many systems. We 
specify neither the physical nature of the oscillators, nor the cou- 
plings necessary for a resonance to occur. In the context of kHz 
QPOs, the oscillators could be taken to be eigenmodes of the ac- 
cretion disk or of an associated structure, perhaps ones already 
discussed in the extensive theoretical literature cited above. We 
do not address the crucial question of the excitation and damp- 
ing of the modes in the turbulent disk, particularly when MRI- 
induced turbulence is present (in the framew ork of two couple d 
oscillators this issue was adressed already bv lVio et~aT1 d2006l) ). 
Eventually, one would hope to identify such resonances in nu- 
merical models of accretion flow. It is only recently that in ertial 
modes have started showing up in numerical simulations dKato 



2004bt lBrandenburg|2005t lArras et al .120061: iRevnolds & Milled 
2008t iMachida & Matsumotol 120081) . It" may be premature to 
look for a numerical validation of the idea of resonance in accre- 
tion disks, at least until two distinct frequencies in a 3:2 ratio are 
identified with some regularity in the simulations. Nevertheless, 
numerical work may already be offering new insights. 

In what follows we review some general properties of inter- 
nal resonance and investigate the amplitude evolution of the os- 
cillations by modeling a slow passage of the frequencies through 
the 3:2 resonance. The results presented in this paper are more 
relevant to neutron-star QPOs where the frequencies wander 
over hundreds of hertz on the time-scale of several hours, than 
to black hole QPOs, whose frequencies are more stable. 



2. Mathematics of the internal resonance 

We model the QPOs using two coupled nonlinear oscillators 
governed by the equations 



x ( + ar c x t = L0 2 c ft{xc, x u , x e , x„), 
X u + <ji- u X u = L0 2 J u {x t , x„, x e , x u ). 



(1) 

(2) 



The oscillations underlying the upper and lower QPOs are iden- 
tified with the dimensionless solutions x u (t) and X((t). The cou- 
pling functions ft and f u are treated as a perturbation, and are 
further assumed to contain only nonlinear terms and to be in- 
variant under time inversion (e.g., their Taylor expansions start 
with the second order and do not contain odd powers of the time 
derivatives if, i„). The coefficients oj( > 0, and <±> u > are the 
eigenfrequencies of the two oscillators. 

In order to study the most relevant resonant effects, we 
concentrate on the case when u> u /u>e is close to 3:2. We 
chose this resonance as it i s most prominent in the QPO data 
dAbramowicz et al.l l2003at iTorok et al.l |2008l) . The degree of 
closeness to the rational ratio is expressed by the detuning pa- 
rameter cr = 2(l>„ - 3aj(. 

Our analysis is applicable when the amplitudes of oscilla- 
tions and the detuning parameter are small and are such that 
ac ~ a u ~ e «: 1 and cr ~ e 2 respectively. In that case, per- 
turbation techniques can be used to find approximate solutions 
of the equations ([TJ and ©. Using the method of multiple scales 
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Fig. 1. Regions where a resonance is present (resonance 
'tongues') in the plane of the eigenfrequency ratio and the to- 
tal energy of oscillations. Each point in the plane corresponds to 
one single {y, £)-disk and different regions to different topologies 
of these disks (see Fig. [2j. Regions (A), (B) and (C) correspond 
to zero, one, and two fixed points in the (y, £)-disk. The parame- 
ters of the system are given in the text. 



(e.g. iNavfeh & Mool 19791) we find the lowest order real dis- 
placements in the form: 



x ( (t) = % [A t (t)e^ ! ] , x u (t) = % [A u (t)e iw "'] . 



(3) 



The higher-order terms add only higher harmonics to the 
quasiperiodic signal described by these solutions. The main dif- 
ference with respect to the solutions of a linear system are that 
the complex amplitudes 



A((t) = X -a t (t)e l ^'\ A„(t) = l -a u {t)e iM) 



(4) 



depend on time. We will refer to the actual frequencies of the 
solution, a>* and a>*, as the "observed frequencies." They are 
shifted with respect to the eigenfrequencies by the corrections 
Aa>e — 4>c and Aa> u — <p u . 

We assume energy conservation, so the two real amplitudes 
may be parameterized using the fractional energy £(f): 



a 2 t =£E, a 2 u = (l-?)-, 0<£<1, 



(5) 



where E = a 2 , + va 2 is the total dimensionless energy and 
v is a constant that d epends on the properties of the system 
dHorak & Karasl 120061) . Similarly, we introduce the phase func- 
tion y(t) as 



7 = 2 <f>u ~ 30f + a-t. 



(6) 



The difference 2a>* - 3co* expressed using the frequency correc- 
tions is 



2u* - 3u* — cr - 34>i + 20„ = y. 



(7) 



Hence, the observed frequencies are commensurable whenever 
y — 0, even if the eigenfrequencies are not exactly so (the devi- 
ation is proportional to e 2 ). 
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The dynamics of the system is then governed by the two first- 
order ordinary differential equations 



1 



^=-^„(l-^£) 3/2 sin r 



1 

7 = cr + - OJ u E 



HS+—0--& + -7I3&- 50(m m cos 7 
v 4 



(8) 
,(9) 



where ft, fi u and fi[ are as yet unspecified dimensionless con- 
stants characterizing a given system. 

The values of the system parameters follow directly from 
the coefficients of the Taylor expansions of the functions and 
f u and depend on the particular physical model. For example, 
iRebuscol d2004l) studied small nonlinear epicyclic oscillations in 
nearly geodesic orbits of test particles in Schw arzschild space- 
time u sing the method of multiple scales. Later. IHorak & Karasl 
d2006l) generalized her approach to motion in an arbitrary ax- 
isymmetric gravitational field. In both works the radial and 
vertical oscillations of the particles are governed by equations 
of the form of CP an d ©, e.g., equations (14) and (15) in 
IHorak & Karasl d2006l) . Moreover, the latter work contains ex- 
plicit formulae for the constants /?, fi u and fif in terms of deriva- 
tives of the effective potential [see the relations after their equa- 
tion (68)]. In addition, a general form of the equations (Q~|) and 
(0 is suitable for s tudying reson ance phenomena in different 
physical conditions. IHorakI d2008l) examined nonlinear interac- 
tion between two epicyclic modes of a slender accretion torus; 
sequation (6) in his paper represents a particular example of the 
general equations discussed here. 

Both, <f(f) and y(t), are slowly varying functions of time, they 
describe a slow modulation of amplitudes and phases during the 
oscillations. The amplitude modulation occurs on the character- 
istic time scale f mo d = e~ 3 f sc- with f osc = being the char- 
acteristic timescale of the oscillations. The two amplitudes are 
anticorrelated and the variations of £ correspond to an exchange 
of energy between the modes that keeps the total energy of the 
oscillations, E, constant. Since the oscillations are nonlinear, the 
observed frequencies are periodically modulated reflecting the 
modulation of the amplitudes. In th e context of QPOs this pro- 
cess has already been discussed by IHorak et"ai1 d2004l) in con- 
nection to the correlations between the QPO fre quencies and the 
phase of the normal-branch oscillations (NBOs) IYu et al.l d200ll) 

The time-evolution of the system can be studied directly 
from equations ([8]) and (0 or with the aid of another integral 
of motion 



*Z + lilt + —(2 - £) + ft! - £){£E) V1 cos 7 
E v 



(10) 



Different solutions trace different curves F(y, £) = const in 
the (y, ^)-disk. For given values of the parameters f3, fi{, fi u , E 
and cr there may exist several fixed points in this disk, each rep- 
resenting strictly periodic oscillations of the system. They cor- 
respond to oscillations whose amplitudes are not modulated and 
whose phases <p((t) and 0„(O are linear functions of time that 
adjust the frequencies of oscillations to the exact 3:2 ratio. The 
fixed points are either extrema or saddles of the functions F(£ , y) 
and their positions in the (7, £)-disk are given by the conditions 
£ = 7 = imposed on equations (O and d9|. Generally, their 
7-coordinate equals kn with k being an integer. 

The number of fixed points in the (7, £)-disk depends on the 
energy E and on the detuning parameter cr. For a fixed value of 
u) u , the latter quantity is given by the ratio of the eigenfrequen- 
cies, % = u> u /a>(. Hence, it is possible to indicate regions in the 
(ft, E) plane for which different number of fixed points exists. 



Figure Q] shows these regions for the system described by the 
fiducial set of parameters,/? = fig = fi u = 1, v = 9/4. 

Some examples of the three possible topologies of the (7, £)- 
disk are given in Figure|2] All of them contain circulating trajec- 
tories for which y(f) monotonically increases or decreases in the 
full range (e.g., < 7 < 2n). The ratio of observed frequencies 
can be expressed (with an accuracy of up to e 3 ) as 



1 + 



2w„ 



(ID 



Hence, this ratio is always greater or smaller than 3/2 if the sys- 
tem evolves along the circulating trajectories with increasing or 
decreasing 7, respectively. 

In addition, closer to the resonance, the (7, £)-disks contain 
a fixed point that corresponds to an extremum of the function 
F(y, £). This point is always surrounded by librating orbits with 
two turning points where 7 changes its sign. When the system 
follows these trajectories the observed frequency ratio oscillates 
around 3:2 . The regions of librating and circulating orbits are 
divided by a separatrix. When the (7, £)-disk contains two fixed 
points, the second point is a saddle point, in which the trajectory 
crosses itself. Consequently, it takes infinite time for the system 
to reach the saddle point from any other point on the separatrix. 



3. Slow passage through the resonance 

Let us now suppose that the eigenfrequencies are tunable and are 
subject to small slow changes. More specifically, we consider the 
case when the eigenfrequencies and couplings are slowly chang- 
ing and the system gradually passes through the resonance. We 
concentrate on that part of this process in which the eigenfre- 
quencies are close to the rational 3:2 ratio. We assume that the 
total relative changes of the eigenfrequencies are small (of the 
order of e 2 ), and the timescale for this change, T, is much longer 
than any other timescale connected to the nonlinear oscillations 
(in particular T » t mo d). 

In general, the oscillations of the system are the solutions to 
equations ([TJ and d2), where a>( and u> u as well as the nonlinear 
functions fc and /„ are known functions that depend explicitly on 
time. Since the characteristic timescale of the process T is much 
longer than t mo[ j, the only change in the equations ©-(O is that 
the eigenfrequencies u>c and aj„ and the coefficients /3, fie, /U„ and 
v become slow functions of time. Close to the resonance, the os- 
cillations are most sensitive to the change of the detuning param- 
eter <x(f) = 2u> u (t) - 3u>((t), while they are practically unaffected 
by small changes of the other quantities. This is because the co- 
efficients of the other terms in equations ([8]l and (0 are always 
multiplied by small quantities (of the order of ~ e 2 and ~ e 3 ). 
We do not consider the "singular" case when the two eigenfre- 
quencies are nearly correlated as oj u as 1.5wf. Because the char- 
acteristic timescale of passage through resonance is taken to be 
much longer than that of the resonant modulation we may ne- 
glect time derivatives of cr and E in the evolution equations (8) 
and (9) that describe the evolution of the system on the timescale 

t mod* 

As the system passes through the resonance, the topol- 
ogy of the (7, £)-disk gradually changes. The system trajectory 
[7(f), £(t)] does not strictly follow the curves F(y, £) = const be- 
cause the ratio cr/E ~ 1 is changing significantly. Indeed both 
quantities are of the order of e 2 and the change of the total energy 
E is negligible because the relative change of v is proportional 
to e 2 . Accordingly, in the following, we take v = const. 
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Fig. 2. Three possible types of topologies of the (y, ^)-disk denoted in Fig.Q]as A, B and C (left, middle and right, respectively). The 
polar angle and the distance from the origin represent the phase function y and the fractional energy respectively. The values of 
the total energy and detuning parameter corresponding to these examples are (cr/u) u , E) = (0.005, 0.1), (0.025, 0.1) and (0.02, 0.1). 
Each curve corresponds to a solution of the governing equations dS) and (0. The disks corresponding to (cr, E) from the resonance 
tongues contain both the librating (white regions), and circulating (shaded regions) trajectories. The fixed points P ± correspond to a 
saddle point or extrema of the integral of motion F(^, y). 




Fig. 3. Behavior of the resonant oscillations of the system undergoing a slow secular change of its eigenfrequencies. The detuning 
parameter <x increases linearly with time so that the system gradually passes through the resonance. The left panel shows the 
trajectory of the system projected into the (y, £)-disk. The right panel shows behavior of the amplitude difference, 6a = a c - a„, 
(measured in arbitrary units) as a function of the ratio the observed frequencies, ft* = to* fait . Because of the slow changes in the 
eigenfrequencies, the energy exchange between modes is not exactly symmetric and the total energy is redistributed between the 
modes when the system passed the resonance. 



The general effect of the resonance transition is illustrated 
by the example shown in Figure [3] The eigenfrequencies are 
changed with constant rate, <x(f) = <n+&t. The left panel shows 
the evolution of the oscillations in the (y, £)-disk and the right 
panel shows the observable quantities (i.e., the ratio of observed 
frequencies ft* versus the amplitude difference da = ac - a u ). 
These Figures are for en = -0.1 w„, & — 3 x 10~ 5 tol, where 
cj u is the upper frequency in the exact resonance. The system is 
described by f3 = —yn = fi u = 5 and v = 9/4, and the initial 
conditions are £ = £(0) = 0.5, y\ = y(Q) = 0. The total energy 
is fixed at E = 0.1. We also tried different rates of passage and 
different initial conditions. 



In all cases the behavior of the system is qualitatively simi- 
lar: the initial value of the detuning parameter corresponds to a 
point (en, E) outside the resonance tongue. Hence, the (y, £)-disk 
initially contains only circulating orbits. One of them is closely 
followed by the system (the outer "circles" in the left panel in 
Figure 3). The ^-coordinate slightly oscillates in correspondence 
to a weak energy exchange between the modes. As the system 
approaches the resonance, the magnitude of these oscillations 
increases gradually. Then the system enters into the resonance 
tongue and the topology of the (y, £)-disk suddenly changes - 
the librating trajectories appear. When the trajectory of the sys- 
tem crosses the separatrix and enters the region of librating tra- 
jectories, the energy exchange becomes asymmetric - i.e., the 
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inflow of the energy into one of the modes exceeds the outflow. 
The system may escape from this region passing the separatrix 
again before it makes a single orbit (as in Fig. 0) or it can be 
trapped in this region for longer time. A trapping in the region 
of librating orbits is accompanied by several changes of the sign 
of y, and therefore by changes of the direction of orbiting in the 
(y, £)-disk. The ratio ft* of the observed frequencies oscillates 
around 1.5 [as follows from equation (fTT~b1 . When the system 
leaves the resonance tongues and as it gets farther away from 
them, the energy exchange between modes becomes more sym- 
metric and less apparent. The final fractional energy £f after the 
passage is generally different from the initial one the energy 
is redistributed due to the resonance. 

We note that the behavior of the amplitudes with respect to 
the frequency ratio as shown in the right pane l of Figure [3] is 
similar to that reported for several atoll sources dTorokKOOai 

An example of a trajectory trapped in the region of the librat- 
ing orbits is shown in Figure [4] The parameters of the solution 
are <x = -0.2w„, & = 4 x lCT 5 ^, p = fi( = -/i u = 5 and 
v = 9/4. The initial conditions are = 0.9, j\ — and the to- 
tal energy is £ = 0.1. Obviously, the nonlinearities adjust the 
observed frequencies so that their ratio oscillates around the res- 
onance value 3:2 during the trapping period. 

We examined the stability of this process against changes 
of the initial conditions and other parameters characterizing the 
passage through the resonance. The effects of changing the ini- 
tial conditions are summarized in the left panel of Figure [3] 
Again, the parameters of the system were set to the values f3 = 
-ju c — ji u — 5 and v = 9/4 and the initial detuning parameter 
and the rate of change were cr o = -0.1 u> u , & = 3 x 10~ 5 «„. The 
solutions trapped in the librating regions and showing for some 
time an episode of what could be termed 'frequency-locking' 
are clearly separated from the others. These solutions are quite 
common for small £ In that case they are accompanied by sig- 
nificant exchanged energy. As the initial fractional energy £ i in- 
creases they become rather exceptional. 

Although different for the solutions with and without 
frequency-locking, the £-£t dependence is fairly independent of 
the initial phase function j\ as it can be seen from the left panel 
of Figure 5, where the solutions with different j\ correspond to 
different symbols. Similarly, we also checked that the rate of the 
transition & has no significant influence on the resulting £ f 
plot. 

The right panel of Figure [5] shows the effect of changing the 
value of the parameter/?. The redistribution of energy in the solu- 
tions without the 'frequency-locking' period becomes more pro- 
nounced as the value of fi increases. In addition, the solution 
with the 'frequency locking' period become more frequent (see 
Table [T). The value /? = corresponds to the case when there is 
no resonant coupling between the two oscillations, and £ = £f . 
Of course, there are no 'frequency-locked' solutions in this case. 

The dependence on the remaining parameters could be dis- 
cussed in a similar way; however, it can be shown that for suf- 
ficiently slow transitions of the resonance the whole dynamics 
depends only on a single parameter, /? = fi^Ej{[i{ - /j.„/v). 
Therefore one does not expect any qualitatively new behavior 
as the other parameters are changed. 

4. Discussion and conclusions 

If the twin kHz QPOs in neutron-stars correspond to two physi- 
cal oscillators that have a nonlinear coupling, the two QPOs may 
exhibit resonance phenomena when the QPO frequencies are 



Table 1. Occurrence of solution with frequency-locking period 
in sets of 200 solutions for different values of j\ and /?. In each 
set the initial value £ is uniformly distributed between and 1, 
while the other parameters are kept constant. 







fi = 


fi = 3 


fi = S 


13= 10 




= 





44 (22%) 


64 (32%) 


89 (48%) 


n = 


= n/2 





47 (24%) 


60 (30%) 


95 (48%) 


7i 


= 71 





43 (22%) 


60 (30%) 


88 (44%) 



close to a ratio of small integer numbers. Among the most char- 
acteristic observable imprints of the internal resonance are (i) 
Exchange of energy between the oscillations (ii) equality of the 
quality factors of oscillations and (iii) locking of the frequency 
ratio of oscillations. Our work was devoted to the mathematics 
of the first and third phenomena. 

Some models of QPOs, which may support an internal res- 
onance of the type described in this work have been described 
in Section 1. In this work we explored only the 3:2 resonance 
which is the one most pro minent in the data dAbramowicz et all 
12003 at ITorok et ai1l2008l) . The examination of other resonances 
may be done in an analogous way. 

Assuming a slow (secular) change of the eigenfrequencies of 
the two coupled oscillators, we demonstrated that the energy is 
redistributed between the two oscillators when the system passes 
through a resonance. This effect can be seen in the change of the 
sign in the amplitude difference when the ratio of the eigenfre- 
quencies ratio crosses a rational value, such as 3:2. This result is 
in qu alitative agre ement with the observations for some neutron 
stars dTorokll2008l) . 

A class of solutions exhibits a time interval during which 
a frequency-locking of the two oscillations occurs. This corre- 
sponds to small oscillations of the actual frequencies around 
the ratio 3:2 over a (greater) range of eigenfrequency ratios 
(Figure 4). 
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Fig. 4. Example of a trajectory with a frequency-locking episode. Left: Trajectory projected onto the (y, £)-disk. frequency locking 
corresponds to the redistribution of the energy when the system is trapped in the region of the librating orbits. Right: The frequency- 
locking effect in the plot showing ratio of the observed frequencies versus eigenfrequency ratio. 




Fig. 5. Energy exchange during the forward transition through the 3:2 resonance (e.g., from the lower to the higher eigenfrequency 
ratios). £ and £f are the fractional energy in the lower oscillation mode before and after the transition. Left: Dependence on the 
initial condition (different symbols denote different initial phases, j\ = 0,n/2 and n) for ft — 5. Right: Dependence on the value 
of the parameter fi (the initial phase is j\ = 0). The branches corresponding to solutions with (denoted by 'FL') and without the 
frequency-locking period are clearly separated. In this plot, £ is sampled uniformly — an indication of the probability of finding the 
system in a 'frequency-locked' solution is given by the relative number of points on the FL branch. 
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